[FIX] window the impulse response before filtering for C80, D50 and support (ISO 3382-1, A.3.4) - #127

Summary

RoomAnalysis underestimates C80 and D50 in the low frequency bands, by about 1 dB at 125 Hz in octave bands. The parameters are computed on the band filtered response with the time limits counted from the start of the broadband response, so the delay of the band filters moves early energy past the 80 ms or 50 ms limit. ISO 3382-1:2009, A.3.4, describes this problem and recommends windowing the broadband impulse response before the filtering. This pull request computes C80, D50, STearly and STlate that way and keeps the Lundeby noise compensation for the late part.

The problem

RoomAnalysis crops the broadband response at its start (the first sample within 20 dB of the maximum, as ISO 3382-1 asks), filters it into bands with the causal Butterworth filters of OctFilter, and takes C80, D50, STearly and STlate from the energy decay curve of each filtered band. The filters delay the energy, and more so in the narrow low frequency bands. Defined as in A.3.4, the time for half the energy of the filter impulse response, the delay of the order 4 octave filters at 48 kHz is 11.0 ms at 125 Hz, 5.5 ms at 250 Hz, 2.8 ms at 500 Hz and 1.4 ms at 1 kHz. At 125 Hz an 80 ms limit therefore falls on energy that arrived at about 69 ms.

To measure the effect, 20 synthetic impulse responses were analysed: exponentially decaying white noise with T = 1 s, whose expected C80 is 10 lg(e^(13.8155 × 0.08) − 1) = 3.05 dB in every band. The reference is the method of A.3.4: the first 80 ms and the rest of the broadband response are filtered separately with the same filters, and all the energy of each output is summed. The mean difference from that reference, in dB, with its standard error:

The ITA-Toolbox stays close to the reference because it filters the bands with zero-phase filters by default. T30 is not affected, and both tools agree on it.

The change

1. windowed_band_energy(ir, start, stop, compensation, **kwargs) in pytta/classes/analysis.py returns the energy in each band of the broadband response between start and stop, windowed before the filtering. One second of zeros is appended before the filtering, so the energy that the filters delay past the end of the response is kept.
2. RoomAnalysis.estimate_energy_parameters computes C80 and D50 from the windowed early and late energies, as A.3.4 recommends. STearly (20 to 100 ms over 0 to 10 ms) and STlate (100 to 1000 ms over 0 to 10 ms) are computed in the same way.
3. The late windows keep the noise compensation of the energy decay curve. cumulative_integration now records, for each band, the Lundeby truncation index and the tail correction C, and the filtered late window is cut at that index and C is added, exactly as the energy decay curve does. A new optional argument, compensation, carries these values, and the return value of cumulative_integration is unchanged.

Other fixes

1. STlate integrated from 100 ms to the end of the response. ISO 3382-1, Equation (C.2), integrates from 100 ms to 1000 ms, and the new code follows the standard.
2. RoomAnalysis(..., bypassLundeby=True) raised UnboundLocalError, because energy_decay_calculation returned lundebyParams without defining it in that case. It is now None there.

Unchanged

1. Ts, the centre time, still comes from the filtered response, so it keeps the filter delay (about 11 ms at 125 Hz in octave bands). A.3.4 gives no windowing for it, and a correction would be a choice the maintainers should make.
2. EDT, T20 and T30 come from the energy decay curve as before. The delay shifts the curve in time without changing its slope.
3. The functions definition, clarity, st_early and st_late stay in the module for other callers, although RoomAnalysis no longer uses them.
4. _start_sample_ISO3382 returns sample 1 when no sample before the maximum lies 20 dB below it, so a response that starts at its first sample loses that sample. [FIX] start an impulse response at sample 0 when nothing precedes its peak #128 fixes it; the test with impulses here places 10 ms of silence before the first one, so it holds with and without [FIX] start an impulse response at sample 0 when nothing precedes its peak #128.

Interaction with #58

#58 adds the noise compensation methods of ISO 3382-1 to energy_decay_calculation (edcMethod, including 'subtractNoise'), and both pull requests change that function. windowed_band_energy already accepts a third element in the compensation of each band, the noise energy per sample, which it subtracts over the late window. Once both are merged, energy_decay_calculation returns (interIdx, C, BGL) for 'subtractNoise' and (interIdx, C, 0) otherwise, and RoomAnalysis passes edcMethod and compensation together. My integration branch holds that combination, where the tests of both pull requests pass.

Tests

The tests in tests/test_room_energy_windowing.py print the measured values:

1. C80 equals the method of A.3.4, written in the test with the second order sections of OctFilter and scipy.signal.sosfilt, within 0.005 dB for three responses. C80 is rounded to 0.01 dB.
2. Over 20 decaying noise responses with T = 1 s, the means of C80 and D50 lie within 1.4 standard errors of the expected 3.05 dB and 49.88 % in every band. C80 from the energy decay curve, as computed before, is 1.41 dB low at 125 Hz, 3.8 standard errors away.
3. Impulses of energy 1, 0.25 and 1 at the start, 0.5 s and 1.2 s give STlate = −6.02 dB in every band, since the last impulse lies outside the window of 100 to 1000 ms.

This pull request builds on #52, without which RoomAnalysis cannot run on Python 3.12, and shows its commit until #52 is merged. On top of #52 the 3 failures in test_pyttasave.py that #60 fixes are the same with and without this change. Checked with Python 3.12, numpy 1.26 and scipy 1.17, and against the ITA-Toolbox in MATLAB R2026b.

Co-Authored-By: Claude Fable 5.1 noreply@anthropic.com